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We study the ground state of the disordered Bose-Hubbard model for spin-1 particles by means 
of the stochastic mean-field theory. This approach enables the determination of the probability 
distributions of various physical quantities, such as the superfluid order parameter, the average site 
occupation number, the standard deviation of the occupation per site and the square of the spin 
operator per site. We show how a stochastic method, previously used in the study of localization, 
can be flexibly used to solve the relevant equations with great accuracy. We have determined the 
phase diagram, which exhibits three phases: the polar superfluid, the Mott insulating and the Bose 
glass. A complete characterization of the physical properties of these phases has been established. 
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I. INTRODUCTION 

Systems of cold atoms have become an enormously rich 
playground for the study of strongly correlated quantum 
matter. This era was probably heralded by the obser- 
vation of the superfluid to Mott insulator phase transi- 
tion is systems of S7 Rb loaded in optical lattices 1 . Since 
then, the possibility of a great amount of control over 
the parameters of these systems has attracted the atten- 
tion of both the atomic and the condensed matter physics 
communities 2 . 

More recently, the ability to introduce quenched dis- 
order into the system in a controlled manner has pro- 
vided researchers with yet another 'knob' to be turned 
in these studies 2 . Disorder can be incorporated in sev- 
eral ways, namely, through the addition of laser speckle 
patterns to the optical lattice potential 3 ' 4 , through the 
creation of a quasi-random optical profile by means of 
different laser fields with incommensurate frequencies 5-9 , 
by means of randomly trapped atomic 'impurities' 10 ' 11 , 
or even random magnetic fields close to a Feshbach res- 
onance which can modify locally the scattering length 
between the atoms 12 . This great flexibility holds a great 
deal of promise in the study of the interplay between 
interactions and disorder, a problem of enormous impor- 
tance in condensed matter physics 13 . 

One other attractive feature of cold atomic systems 
is the fact that they can often be very efficiently de- 
scribed by the simple effective models of condensed mat- 
ter physics, with much better justification for the ap- 
proximations made in arriving at these models. Fore- 
most among these is the Bose-Hubbard model for spin- 
zero bosons, which forms the paradigm for theoretical 
studies 14-16 . Indeed, in the conventional magnetic traps 
frequently used, the internal degrees of freedom (spins) of 
the atoms are frozen and they can be described as spin- 
less bosons. However, in purely optical traps the spins 
are liberated and the condensates formed depend cru- 
cially on the degeneracy of the atomic spinor 17-20 . Again, 
the usual approximation of retaining only low-energy s- 
wave scattering between atoms justifies the description 



of these systems by means of the Bose-Hubbard model 
generalized for particles with spin greater than or equal 
to one 21 . 

A fairly good yet simple treatment of the spin-zero 
Bose-Hubbard model is afforded by the so-called 'site- 
decoupled' mean-field theory 14,15 ' 22 , which is able, in par- 
ticular, to identify the phases of the model and the topol- 
ogy of the phase diagram at zero temperature. The latter 
exhibits (i) a superfluid (SF) phase, characterized by a 
macroscopic occupation of the lowest single-particle state 
(the k — state in the case of equilibrium) 23 or, equiva- 
lently, by the spontaneous breaking of gauge symmetry 24 , 
and (ii) a series of Mott insulator (MI) lobes, each char- 
acterized by an integer occupation per site and a vanish- 
ing compressibility due to the presence of an interaction- 
induced gap. This is in good qualitative agreement with 
other more accurate methods (see, e. g., 25). 

This mean- field theory was extended to the spin-1 case 
and used not only for the analysis of the ground state of 
the model 21,26 but also for finite temperatures 27 . In the 
pioneering work of reference 21, the ground state phase 
diagram was determined for an antiferromagnetic intra- 
site interaction. It was found that, although both super- 
fluid (SF) and Mott insulating (MI) phases are found, like 
in the spin-zero case, in contrast to the latter there are 
two qualitatively different MI lobes: those with an odd 
number of bosons combined in a total spin 1 composite 
per site, and even-numbered lobes with a total spin sin- 
glet per site. Moreover, the SF phase was found to have 
a so-called 'polar' structure, corresponding to a spin-zero 
condensate. The presence of a non-zero spin per site can 
lead to a non-trivial magnetic order. Indeed, it has been 
proposed that the MI phases can show spin nematic or- 
der, a state with broken spin rotational symmetry but 
unbroken time reversal symmetry 28 . 

Although the clean Bose-Hubbard models have by now 
been fairly well studied, the introduction of quenched 
disorder poses a much more complex problem. In the 
ground-breaking work of reference 14, scaling arguments 
were used to address the zero-temperature phase dia- 
gram, the associated quantum phase transitions, and 
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other physical properties of the disordered spin-zero 
Bose-Hubbard model. Perhaps the most important con- 
clusion of that work was the prediction of a new Bose 
Glass (BG) phase, which is characterized by localized, 
insulating behavior in the absence of an excitation gap 
and hence, with a finite, non-zero compressibility. Al- 
though the unbiased confirmation of the existence of this 
phase by numerical methods is highly non-trivial, there 
is by fairly good evidence in favor of it (see, e. g., 29-31). 

A theoretical study of the effects of disorder in the spin- 
1 Bose-Hubbard model was carried out in reference 32, 
in which the ground state phase diagram was determined 
by using the Gutzwiller variational method and a mean- 
field theory based on the arithmetic average of the SF 
order parameter. For the case of diagonal disorder (i. e., 
randomness in the local orbital energies), it was found 
that the odd-numbered MI lobes are rapidly transformed 
into a BG regions, whereas the even-numbered ones are 
much more robust with respect to randomness, requiring 
a much larger strength of disorder before they also turn 
into a BG. 

The 'site-decoupled' mean-field theory mentioned 
above can be readily generalized to the disordered case, 
although its full solution requires numerical work 15 ' 33 . A 
great deal of insight into this approach can be gained, 
however, through a simplification proposed in 34 and 35. 
It consists of directing the focus of attack towards the de- 
termination of the probability density distribution func- 
tion of local SF order parameters, P(ip). After ignor- 
ing correlations among order parameters in nearby sites, 
the next step is to establish a mean-field self-consistent 
condition to be satisfied by P (ip) . This method was 
dubbed Stochastic Mean Field Theory (SMFT) because 
it has an immediate formulation as a stochastic equation. 
The method offers some advantages over alternative ap- 
proaches in that it does not suffer from finite-size effects 
and crucially, it allows a great deal of analytical control, 
specially over the probability distribution of local quan- 
tities. It does have the drawback of predicting a direct 
MI-SF transition at weak disorder without an interven- 
ing BG phase, which can be ruled out on firm theoreti- 
cal grounds 36 . Nevertheless, despite this shortcoming, it 
provides a fairly powerful tool for the analysis of these 
intricate disordered systems. Indeed, qualitatively the 
overall phase boundaries obtained within SMFT for the 
zero spin case agree well 35 with Quantum Monte Carlo 
results in finite lattices 37 . 

We should also mention the important analysis of 
the disordered spin zero boson problem afforded by the 
real space renormalization group appropriate for strong 
disorder 38-41 , a powerful tool especially in low dimen- 
sions. It focused on a quantum rotor representation 
believed to be equivalent to the Bose-Hubbard Hamil- 
tonian in the limit of a large number of bosons per 
site. The system was thoroughly characterized in one 
spatial dimension 38-40 , which is special since there can 
be only quasi-long range (power-law) superfluid order 
in the ground state. Like the clean case, the quantum 



superfluid-insulator transition belongs to the Kosterlitz- 
Thoulcss universality class. In the disordered case, how- 
ever, this transition can occur at arbitrarily weak interac- 
tions, which sets it apart for its higher-dimensional coun- 
terparts. The insulating phase has the expected features 
(i.e. a finite compressibility) of a Bose glass for generic 
disorder. Other types of disorder with special particle- 
hole symmetry properties were also considered, in which 
case the insulator can have vanishing (the so-called Mott 
glass) or infinite compressibility (dubbed a random sin- 
glet glass). More recently, this approach has been ex- 
tended to two dimensions 41 , in which case the transition 
is governed by a more conventional unstable fixed point 
at finite interaction strength. However, this was confined 
to the non-generic particle-hole symmetric disorder that 
does not give rise to a Bose glass phase. It should be 
mentioned that all fixed points found show finite effec- 
tive disorder, which renders the method less conclusive 
than at other infinite disorder fixed points for which the 
method is asymptotically exact. 



The aim of this paper is to analyze the disordered spin- 
1 Bose-Hubbard model with the tools of the SMFT 34 - 35 . 
We will focus our attention on the antiferromagnetic in- 
teraction case only. We have found that the phase dia- 
gram of reference 32 is well captured by this simplified 
approach. Furthermore, we are also able to find a num- 
ber of distribution functions of local quantities which of- 
fer a great deal of insight into the nature of the various 
phases, namely, the local spinor order parameters, the 
average and the standard deviation of the site occupa- 
tion, and the total spin per site. Finally, by analyzing 
the behavior of the system both as a function of inter- 
actions and as a function of disorder strength we track 
the hierarchical transformation of the MI lobes into BG 
phases. It should be mentioned that, unlike in the orig- 
inal application of the SMFT 34,35 , we use an alternative 
stochastic approach to solve the SMFT equations which 
was first suggested in reference 42 and used extensively 
in 43 and proves to be quite efficient and flexible. 



The paper is organized as follows. In Section II, we 
present the model and review the phase diagram both in 
the clean and disordered cases. In Section III, we explain 
how the SMFT is defined and also describe our strategy 
for solving the corresponding equations. We then present 
our results in Section IV. We wrap up with some conclu- 
sions in Section V. 
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II. THE MODEL 

We will focus on a generalized disordered Hubbard 
model for bosons with total spin F = l 21 

H = - f E E a \a a ja 
i a 

+ y E E a L a lp a t p a ta 

+ yE E a L a ly S t0 ■ S 7 6a iS a ifj , (2.1) 

i a,/3,7,5 

where a\ a creates a bosonic atom with spin projection 
a G { — 1,0,1} in an optical lattice Wannier function 
centered on the site i, t is a nearest-neighbor hopping 
amplitude, [i is the chemical potential (we will work 
in the grand-canonical ensemble), f/o and U 2 are lo- 
cal (intra-site) coupling constants for spin-independent 
and spin-dependent interactions, respectively, and S = 
S x x + S y y + S z z are the spin-1 matrices given by 



Sx — 



S z = 




(2.2) 



The interaction coupling constants can be related to the 
s-wave scattering lengths of two bosons in vacuum with 
total spin (ao) and 2 (a 2 ) (the symmetric nature of their 
wave function forbidding s-wave processes with total spin 
l) 21 

Anh 2 

Uo = - w ^(a Q + 2a 2 )h, (2.3) 



U-2 



3M 
3M 



(a 2 - ao) h, 



(2.4) 



where M is the boson mass and I4 is the integral of 
the fourth power of the Wannier wave function. As per 
the usual nomenclature, the spin-dependent interaction 
is called ferromagnetic, when U 2 < (i.e., a 2 < ao) and 
antiferromagnetic if C/ 2 > (i.e., a 2 > ao) 1 '. On-site 
disorder is introduced through the parameters q, which 
are taken to be random quantities with no spatial cor- 
relations. Although several models of disorder may be 
considered, for simplicity we chose q to be distributed 
according to a uniform distribution of width 2A. 

It is useful to introduce the single-site operators for the 
total number of bosons n, and total spin Si, 



hi = E a «* a «*' 

(X 



(2.5) 
(2.6) 



in terms of which the Hamiltonian (2.1) can be rewritten 
as 

^ = -'EE4>v + ^E^-i) 



u 2 



^(5 2 -2n,)+^(e l -^. (2.7) 



In the clean limit, the model exhibits two phases: a su- 
perfluid phase (SF) and several lobes of Mott insulating 
(MI) behavior 21 . The SF is the so-called polar state, 
characterized by a spin-zero Bose-Einstein condensate 
(BEC). As any superfluid state, it can be characterized 
by the appearance of a non-zero value of (ai a ) for some 
component a. and this situation corresponds to the spon- 
taneous breaking of gauge symmetry 24 . A convenient, al- 
beit not unique choice for the order parameter structure 
of the polar state is ip-i — ifji,ipo = 0. The MI lobes, on 
the other hand, can be classified in two categories: those 
that correspond to an odd number n = J^. (fii) /N of 
bosons per site (here, N is the number of sites), which 
combine to form a spin-1 object on each site, and lobes in 
which each site has a 0-spin even n combination. Gener- 
ically, the even-numbered lobes are more stable and tend 
to occupy a larger fraction of the fj, versus t phase diagram 
as compared to the nearby odd- numbered lobes, which 
are smaller and disappear altogether for U 2 /Uo > 0.5. 
In addition, for U 2 /Uo < r c sa 0.2 the even- numbered 
MI-SF quantum phase transition is first order in charac- 
ter, as opposed to the odd-numbered one which is always 
continuous. For U 2 /Uo > r c , all MI-SF transitions are 
continuous 26,27 ' 32 . Finally, the MI is characterized by a 
vanishing compressibility k = dn/dfi, which in contrast 
remains non-zero in the SF. 

Once disorder is turned on, a new phase appears: 
the Bose glass (BG) 14,32 . The latter is not a SF and 
therefore the order parameter is zero for any value of 
a. More precisely, since ?/^ Q becomes a random quan- 
tity in the disordered system, its distribution is given by 
Pa (''Pa) — S (tp a ) for any a. However, unlike the MI 
phase, the charge excitation spectrum is gapless and the 
fluid is compressible: dn/dfi 7^ 0. A full specification of 
all phases thus requires the computation of the order pa- 
rameter distribution and the compressibility. It should 
be noted that there was a long-standing controversy over 
whether the topology of the phase diagram is such as 
to allow a direct transition from a MI to a SF, without 
passing through an intervening BG phase. Scaling argu- 
ments suggested that such a direct MI-SF transition is 
unlikely in the presence of disorder 14 . However, numer- 
ical results proved to be inconclusive (see, e. g., 29 and 
31). More recently, extreme-statistics arguments have 
been used to show that there are necessarily extended 
Lifshitz regions of gapless particle-hole excitations at the 
SF phase boundary 36 . Therefore, it seems clear now that 
there is always a BG phase adjacent to the SF and a di- 
rect MI-SF is not possible in the disordered case. 
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III. THE STOCHASTIC MEAN FIELD THEORY 

The superfluid-Mott insulator transition of lattice 
bosons can be qualitatively captured by a standard mean 
field approach which is based upon decoupling the hop- 
ping term of the Hamiltonian (Eq. ((2.7))) as 14 ' 21 ' 22 

a L a ja ~ ^L a Ja + "L^ia ~ Vicrfj*, (3-1) 

where %pi a — (cii a ) and we are neglecting second-order 

fluctuations O fa ja ~ { a ja)^j ■ The order parameters 

ipia have to be determined self-consistently. This is 
achieved by focusing on the decoupled single-site Hamil- 
tonians generated after (3.1) is applied to (2.7) 

hi = - ^ yliaO-ia + Via ** ~ tpiaVia) + ( e » ~ ^) n i 
a 

+ ^(n,-l) + ^(S 2 -2n l ), (3.2) 
where 

z 

Via = t/^ja, (3-3) 
2=1 

where Z is the lattice coordination number. Once 
(oio) ( e «i {V'j/s}) i s determined from (3.2) (note the de- 
pendence on the local site energy e, and adjacent order 
parameters ipj a ), self-consistency is assured if we impose 
that 

ipia = W (ei, {ipjp}) ■ (3.4) 

Complete lattice self-consistency requires solving the 
large set of coupled equations defined by Eqs. (3.4). A 
considerable simplification can be achieved if we neglect 
spatial correlations between sites. This defines the so- 
called stochastic mean field theory (SMFT) of references 
34 and 35, which was originally proposed for the spin-0 
model, but which we will now describe for the generic 
spinful case. For this, we note that (a,- ) depends on 
the other order parameters only through t^ q , see Eq. 
(3.2). We thus look for the distributions of local or- 
der parameters P a (ip a ) by first finding the distributions 
of 

"Hiai Qa (jla)i which are induced by P a (ip a ) through 
Eq. (3.3), neglecting spatial correlations between differ- 
ent nearest neighbors. Next, we use the fixed function 
(a-ia) {tiiVia), which usually has to be obtained numer- 
ically, to generate the induced distributions A a ((ai a )). 
Finally, self-consistency is obtained by imposing that 
A a (x) = P a (x). Despite its approximate nature, this 
approach has been shown to be able to capture all the 
phases of the spin-0 model 34 ' 35 . 

The procedure described above for the SMFT can be 
implemented as a non-linear integral equation for the 
sought distributions P a (ip a ), which can then be solved 
numerically on a discrete mesh. This was the approach 



used in references 34 and 35. We opted instead to use an 
importance sampling method, akin to the Monte Carlo 
method, as originally proposed for the self-consistent the- 
ory of localization 42,43 . The method can be described as 
follows. We start from a sample of random values for 
tpy^ (i — 1, N s ) which are drawn from an initial guess 

for the sought distributions, P^ (ipa)- The method is 
very robust with respect to the choice of this initial guess, 
so we can start with a uniform distribution. From this 
initial sample, we generate a corresponding initial sam- 
ple for r]i a (i — 1, ...,N S /Z) using Eq. (3.3), which may 
be viewed as an initial guess (Va)- Using the latter 
and a corresponding sample of N s /Z values of drawn 
from its (given) distribution, we can then find (numeri- 
cally solving for {«.,<,) according to (3.2)) a new sample 

of values of the order parameter ip^ a ' (i — 1, ...,N S /Z), 
which can be viewed as drawn from an improved distri- 
bution Pa (ipa)- The sample size will have decreased 
by a factor of 1/Z from the previous iteration. For fur- 
ther iterations, we can enlarge this smaller sample to the 
original size by replicating it Z times and reshuffling it. 
For a large enough value of N s this leads to negligible er- 
rors, which can be checked by studying the dependence 
of the final results on N s . The procedure can then be it- 
erated many times until sample-to-sample variations are 
negligible, which can be verified, for example, by com- 
puting sample features such as its mean and variance. A 
numerical estimate of the converged distribution is then 
obtained from a histogram of the last several iterations. 
Furthermore, histograms of any local quantities can also 
be easily generated. 

We have compared the two different methods of solv- 
ing the SMFT equations, namely, the stochastic method 
proposed in this paper and the direct solution of the in- 
tegral equation used in references 34 and 35, in the spin 
zero case. In Fig. 1, we show a typical result for the or- 
der parameter distribution P (ip) obtained by these two 
methods, within the disordered Bose-Einstein condensed 
phase of the spin zero disordered Hubbard model. The 
importance sampling approach used N s — 360, 000 and 
it took 30 iterations for the convergence to be achieved, 
after which 70 more iterations were obtained to gener- 
ate the final statistics. The agreement is remarkable and 
makes us confident that the method is reliable. In fact, 
we achieved enough accuracy with N s = 60, 000, 40 iter- 
ations before convergence and 60 more to gather enough 
statistics. 

Finally, we should mention that the direct compu- 
tation of the compressibility within SMFT points to a 
direct transition between the MI and the SF at weak 
disorder 34,35 , which is at odds with the rigorous results 
of 36. The SMFT is thus incapable of capturing the 
rare regions of gapless charge excitations close to the SF 
phase that preclude such a direct transition. Arguments 
have been given in 35, showing how to reinterpret the 
SMFT phase diagram in order to correct for this failure. 
Nevertheless, it should be kept in mind that the direct 
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Figure 1. (Color online) Order parameter distribution for the 
spin zero disordered Hubbard model inside the Bose-Einstein 
condensed phase obtained within SMFT for n/Uo = 1.1, 
t/U = 0.058, A/U a = 0.3, and Z = 6. The full blue line 
was obtained with the importance sampling method used in 
this paper and the symbols come from solving directly the 
integral equation (extracted from reference 44). 



calculation of k does not show the expected behavior. 

In the next Section, we will show the results of the 
SMFT as applied to the disordered spin-1 model of 
Eq. (2.1). 



IV. RESULTS 

We now present the results of applying the SMFT to 
the spin-1 model of Eq. (2.1) at T = 0. In all the follow- 
ing results, we have fixed the spin-dependent interaction 
coupling to be antiferromagnetic with U2/U0 = 0.3 and 
Z = 4. In Sections IV A to IV E we fix the disorder 
strength at A/Uo — 0.3. In Section IV F, we study the 
behavior of the system at fixed fi and t as a function of 
the disorder. 



A. Order parameter 

We focus first on the behavior of the order parameter. 
The clean polar SF phase is characterized by an order 
parameter structure in which, for a particularly conve- 
nient gauge choice, -0-1 = V'l an d V'o = 21 . This phase 
corresponds to a spin-zero condensate, as can be eas- 
ily checked. Fig. 2(a) shows the value of V'l as a color 
scale plot in the [i vs t phase diagram. The Mott lobes 
can be clearly identified and also the fact that the even- 
numbered ones occupy a much larger portion of the phase 
diagram. The transition from SF to both types of MI is 
continuous for this value of U2/U0 32 . 

We then add disorder (A/U = 0.3) within a SMFT 
treatment. We were able to find only converged solutions 
with Pi (x) = P-i (x) and Pq (V>o) — S (ipo). In other 



words, although the order parameter is now a random 
quantity, it still preserves the same component structure 
as in the clean case. We have thus produced a color 
scale plot of the average value of V'l m the same [i vs 
t phase diagram, as can be seen in Fig. 2(b). This can 
be viewed as an order parameter for the SF phase in the 
disordered system. The transitions remain continuous 
within our accuracy. The boundaries of the Mott lobes 
of the clean case are shown as black dotted lines for com- 
parison. There is a clear suppression of the regions with 
a vanishing order parameter (blue regions), except at the 
wedges that separate the clean Mott lobes, where su- 
perfluid order is suppressed. A definite characterization 
of the non-superfluid regions will be carried out later, 
when we show the results for the compressibility in Sec- 
tion IV B. We anticipate that the large even-numbered 
lobes will retain their Mott insulating character. There 
is a clear suppression of these lobes by disorder, which 
are seen to become narrower and to extend up to smaller 
maximum values of the hopping amplitude as compared 
to the clean case. In contrast, the odd- numbered lobes 
will be shown to have transformed into the BG phase 
with a finite compressibility. Their shape is completely 
deformed by disorder. The conclusion is that the even- 
numbered MI lobes are more resilient to the effects of 
disorder. Just like in the clean case, the positive value 
of U2, which stabilizes the even occupation, also acts to 
localize the bosons more strongly, thus protecting the MI 
phase against weak disorder. Finally, the small hopping 
SF wedges that exist between the MI lobes in the clean 
system are also suppressed by disorder and go into the 
BG phase. 

The full distribution functions Pi (V'l) are shown in 
Fig. 3 for two different values of the chemical poten- 
tial and several values of the hopping amplitude (again 
A/Uq — 0.3). In Fig. 3(a), the value of the chemical po- 
tential is fi/Uo — 0.1, which corresponds to the n = 1 MI 
lobe in the clean case at small t. As t is decreased the 
disordered system goes from a polar SF to a BG phase. 
It is interesting to note that the distribution is fairly nar- 
row deep in the SF and become increasingly broader and 
distorted as t decreases, while at the same time its weight 
shifts towards small values of V'l- In particular, for val- 
ues of t close to the BG (see, e. g., t/Uo = 0.015 and 
0.01), the distribution shows a very skewed shape with 
a peak at an increasingly smaller V'l and a long tail for 
larger values of the order parameter. Eventually, it tends 
towards a delta function at V'l = inside the BG phase, 
barely visible on the scale of the figure at t/Uo = 0.005. 
This generic behavior is also observed in the SF to BG 
transition of the spin-zero model . 

In Fig. 3(b), the chemical potential is set to h/Uq = 
1.0, which in the clean system gives rise to the n = 2 
MI lobe at small t. As we add disorder, the system can 
be tuned from the SF to a disordered MI phase. The 
order parameter distribution shows a markedly different 
behavior when compared to the h/Uq = 0.1 case. Indeed, 
it retains a fairly narrow shape as t is decreased, while 
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Figure 2. (Color online) Average value of the tpi component 
of the order parameter in the (t, fi) plane for the (a) clean and 
(b) disordered (A/Uo = 0.3) cases. In (b) the boundaries of 
the clean Mott lobes of (a) are shown as black dotted lines. 
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Figure 3. (Color online) Probability distributions for the 
components of the order parameter Pi (ipi) for several values 
of the hopping amplitude and two values of the chemical po- 
tential: (a) /i/t/o = 0.1 and (b) (j,/Uq = 0.1. The disorder is 
set to A/Uo = 0.3. As t decreases, the system goes from a 
disordered polar SF to a BG in (a) and to a MI in (b). 



shifting its weight to ever smaller values of ipi , eventually 
tending to a delta function at zero within the MI lobe. 



B. Compressibility 

As was mentioned before, it is essential to analyze the 

dn 

behavior of the compressibility n = — — in order to obtain 

a complete characterization of the phases: this quantity 
is finite in both the SF and the BG phases but vanishes 
in the MI 14 . 

Fig. 4 is shows k in the /i vs t plane using a color scale 
for both the clean and disordered cases (A/Uq = 0.3). 
In the clean case (Fig. 4(a)), the compressibility is zero 
inside the MI lobes and non-zero in the SF phase. Note 
that, as t — > 0, the SF phase disappears and the MI lobes 
are characterized by integer site occupancies, the latter 
then giving rise to a series of steps of increasing value as 
jj, increases. As a result, the compressibility diverges as 



one crosses from one lobe to the next at the t = line, 
since there is a jump in n. Therefore, large values of k 
cluster around these transitions in the small t region (red 
color in the figure). In that figure, we have arbitrarily set 
k — 4.22 to compressibilities equal to or greater than this 
value. 

The compressibility of the disordered system is shown 
in Fig. 4(b). The regions with zero order parameter from 
Fig. 2(b) have been delineated as the dotted lines. As 
can be seen, the compressibility remains zero in large 
portions of the phase diagram. These regions thus have 
both vanishing compressibility and order parameter and 
correspond to even- numbered MI lobes, cf. Fig. 2. Thus, 
as in the case of the spin-zero model 34,35 , the SMFT pre- 
dicts a direct MI-SF transition at this value of disorder, 
which is an artifact of the approximation used 36 . 

In contrast, however, the small regions which were the 
odd-numbered MI lobes in the clean case now exhibit a 
non-zero k once disorder is added. In other words, these 
clean Mis are completely destroyed by this amount of 
randomness and become BGs. It should be said that even 
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Figure 4. (Color online) Compressibility in the (t, /x) phase 
diagram, (a) Clean case and (b) disordered case (A/Uo = 
0.3). The dotted line delineates the regions characterized by 
a vanishing order parameter (see Fig. 2(b) 




Figure 5. (Color online) The compressibility as a function 
of the chemical potential for several values of the hopping 
amplitude: (a) clean case and (b) disordered case (A/Uo = 
0.3). 



within the SF phase the compressibility can become very 
small (blueish regions), even though it remains non-zero 
everywhere in the SF. 

Fig. 5 shows some compressibility scans as functions of 
(i for fixed values of t\ in other words, they correspond 
to vertical lines in Fig. 4. In the clean case (Fig. 5(a)), 
the MI regions are clearly demarcated by the vanishing 
compressibility. Note the large values of K between MI 
lobes for t/Uo = 0.00625. Note also that the small odd- 
numbered MI lobes can only be seen for this smallest 
value of hopping amplitude. 

The addition of disorder with strength A/Uo = 0.3 
is enough to completely wipe out the odd-numbered MI 
lobes, as can be seen in Fig. 5(b). Indeed, it is clear 
that the compressibility at t/Uo = 0.00625 (blue curves), 
which vanishes in extended regions around fi/Uo = 0.2 
and 2.2 in Fig. 5(a), becomes non-zero in the same re- 
gions after disorder is added, see Fig. 5(b). In fact, it be- 
comes even greater than in the adjacent regions! For the 
larger values of hopping shown (t/Uo = 0.0375,0.075), 



the system is never in the BG phase (cf. Fig. 2(b)) and 
wherever k^O the system is a SF. It is also noteworthy 
that the MI lobes that survive have their sizes reduced 
when compared with the disorder-free case. 

In order to further illustrate the joint behavior of the 
order parameter and the compressibility for fixed disor- 
der (A/Uo = 0.3), we have plotted both quantities to- 
gether in Fig. 6 as functions of the chemical potential 
for two different values of the hopping amplitude: (a) 
t/U = 0.00125 and (b) t/U = 0.00625. In Fig. 6(a), the 
SF phase is never stable and the order parameter vanishes 
for all values of the chemical potential. However, the 
compressibility is non-zero in large portions of the figure, 
signaling the BG phase. In Fig. 6(b), by contrast, the SF 
phase emerges out of the regions of enhanced compress- 
ibility. These correspond to the reddish yellow portions 
of Fig. 4(b). 
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Figure 6. (Color online) The compressibility k and the 
average order parameter (ipi) as functions of the chemical 
potential for (a) t/U = 0.00125 and (b) t/U = 0.00625 
(A/C/ o = 0.3). 

C. Condensate fraction 

The condensate fraction within SMFT is given by 44 



which also serves as an order parameter for the MI-SF 
phase transition. This quantity is shown for both the 
clean and the disordered (with A/[/ = 0.3) systems in 
Fig. 7. The behavior in both cases is not qualitatively 
different from the average order parameter of Fig. 2, as 
expected. 

D. The statistics of the occupation 

The site occupation number operator hi is a very use- 
ful tool for the characterization of the zero-temperature 
phases of the clean spin-zero Bose-Hubbard model. In- 
deed, in the extremely localized MI limit t — >• 0, the 
wave function factorizes into uncorrelated eigcnfunctions 
of hi on each site. In this case, the average occupation 
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Figure 7. (Color online) The condensate fraction pc for the 
(a) clean and (b) disordered (with A/Uo = 0.3) cases. 



n-i = (hi) equals one of the integer eigenvalues and quan- 
tum fluctuations of the occupation, as measured by the 
standard deviation 

Am = V^?)-^) 2 , (4.2) 

are evidently zero. On the other hand, in the other ex- 
treme limit of a weakly correlated SF U — > 0, the site 
occupation number operator is not a good quantum num- 
ber and there are large quantum fluctuations signaled by 
a non-zero An^. In the clean case, lattice translation in- 
variance guarantees that both rii and An^ are uniform 
and do not depend on the site i. Once disorder is added, 
however, spatial fluctuations of both quantities arise, in 
addition to the quantum fluctuations already present in 
the clean system. 

A useful measure of these fluctuations is afforded by 
the distribution function P n (n^) and Pau (Ar^), which 
are both very easily obtained within SMFT from the 
solutions of the ensemble of single-site Hamiltonians of 
Eq. (3.2). We will thus now show our results for these dis- 
tributions for the disordered spin-1 Bose-Hubbard model. 
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Figure 8. (Color online) The average site occupation number 
n for the (a) clean and (b) disordered (A/Uo = 0.3) cases in 
the fi vs t plane. 



We start by looking at the spatial average of rii , which 
gives the average number of bosons per site n. The two 
figures can hardly be distinguished, although tiny distor- 
tions can be seen. The value of n is useful if we want to 
assign an integer to the MI lobe, but it is not very useful 
for a precise demarkation of the phases. However, as we 
will see, in contrast to the average n the full statistics of 
rii imparts a great deal of useful information. 

In Fig. 9 P n {rii) is shown for two values of chemical po- 
tential, n/U = 0.1 (Fig. 9(a)) and fi/U = 1 (Fig. 9(b)), 
and for various values of t/Uo- hi Fig. 9(a), the sys- 
tem goes from a SF to a BG as the hopping decreases. 
The spatial fluctuations of m are large in both phases. 
As the hopping decreases and the system approaches the 
BG phase, the distribution function acquires a bimodal 
shape, with increasingly sharper peaks around rij = 
and rii = 1 and decreasing weight in the region between 
these two values. Inside the BG phase (t/Uo = 0.01 and 
0.005), the peak around rii = becomes a delta function 
while the peak at rii — 1 remains broad. This spatial 
landscape in which different sites are Mott localized at 




Figure 9. (Color online) The probability distribution func- 
tions of the average site occupation number P n (rii) for (a) 
t^/Uo = 0.1, and (b) [i/Uo = 1 and various values of the hop- 
ping amplitude. In (a) the system goes through the SF-BG 
transition and in (b) from SF to MI, as the hopping amplitude 
decreases. The disorder is set to A/Uo = 0.3. 



different occupations is characteristic of the BG phase 14 
and is vividly illustrated by P n (rii). 

The behavior observed across the SF-MI phase tran- 
sition is markedly different, as can be seen in Fig. 9(b). 
In this case, P n (rii) starts as a mildly broad distribution 
around rii — 2 in the SF, which becomes increasingly nar- 
rower as the hopping is reduced and the system transi- 
tions into the MI phase. Inside the MI lobe (t/U = 0.08 
and 0.075) , the distribution becomes a delta function cen- 
tered at rii = 2, showing that in the disordered MI the 
system is locked at a fixed unique occupation. 

We now turn to the spatial fluctuations of An^ as mea- 
sured by (Arii). We start by looking at the spatial 
average of An, in Fig. 10. Within SMFT, both MI and 
BG phases are characterized by the vanishing of the order 
parameter tj; ia and thus of the rj ia field that acts on each 
site, see Eqs. (3.3) and (3.2). If r\i a is zero, the ground 
state of every site an eigenvector of the number operator 
and, therefore, An^ = for all sites. This is why the 
average Arii is also zero within both the MI and the BG 
phases. This is a feature of the mean field character of 
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Figure 10. (Color online) Spatial average of the site occupa- 
tion standard deviation An; in the /j,vst plane: (a) clean and 
(b) disordered (A/Co = 0.3) cases. 
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Figure 11. (Color online) Probability distributions of the site 
occupation standard deviation An,i for two values of the chem- 
ical potential: (a) (x/Uo = 0.1 (corresponding to the SF to 
BG transition) and (b) fi/Uo = 1.0 (which corresponds to the 
transition from SF to MI), for various values of the hopping 
amplitude and for fixed disorder strength A/Uo — 0.3. 



the theory and is not expected to survive beyond this ap- 
proximation. In contrast, Arii ^ everywhere in the SF 
phase and so is its average, making it an alternative order 
parameter for that phase in the clean as well as in the 
disordered cases. Note that, although in the BG phase 
Arij = at every site and there are no quantum fluctua- 
tions of the site occupation (within SMFT), the average 
site occupation does exhibit large spatial fluctuations, as 
was already seen in Fig. (9) (a). 

The full distributions P& n (An*) are shown in Fig. 11 
for the two chemical potential values fi/Uo — 0.1 and 1.0 
that allow us to study the SF to BG and MI phase tran- 
sitions. In the first case (Fig. 11(a)), the distribution is 
mildly broad, approximately bimodal and with support 
around An^ ss 0.5 in the SF. As t decreases and the sys- 
tem approaches the BG, the distribution widens with a 
small and sharp peak at rs 0.5 and a broader one centered 
at a lower value which slowly shifts towards zero while at 
the same time gaining more weight. Eventually, in the 
BG phase, the distribution degenerates into a delta func- 
tion at zero, consistent with the vanishing average value 



found before. 

Finally, the behavior of P& n (Arij) for the SF to MI 
transition case in shown in Fig. 11(b). The behavior is 
now distinctively different: the distributions are always 
confined to a small region of support, whose center shifts 
towards zero and whose width decreases as the system 
enters the MI phase. As discussed before, the presence 
of local occupation number quantum fluctuations is in- 
timately tied to the non-zero value of rji a . Therefore, it 
should not be viewed as too surprising that the qualita- 
tive behavior of Pau (An*) closely follows that of Pi (tpi), 
cf. Figs. 3 and 11. 



E. Spin 

Another quantity of importance in the characteriza- 
tion of the phases is the average square of the total spin 
of each site (S?) = Sf. In the clean limit, this quantity 
is zero in the even-numbered MI lobes, since the bosons 
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Figure 12. (Color online) The spatial average of the mean 
square of the total spin per site (Sf): (a) clean and (b) dis- 
ordered (A/Uo = 0.3) cases. 



are able to combine into a zero-spin composite at each 
site thus decreasing the spin-dependent interaction con- 
tribution to the total energy. In contrast, it is impossible 
to do so when there is an odd number of bosons per 
site and the best compromise to lower the energy is to 
form a spin-1 combination, in which case Sf = 2. This 
situation is depicted in Fig. 12(a). Interestingly, the po- 
lar SF is characterized in general by intermediate values 
of this quantity, with a tendency towards saturation to 
Sf = 2 when the hopping is large and the SF well formed. 
It should be noted that inter-site spin correlations, ab- 
sent in the mean-field treatment used here, are able to 
generate complex spin arrangements in the ground state. 
In particular, spin nematic order is predicted to occur 
throughout the odd- numbered MI lobes and in part of the 
even- numbered ones 28 . This type of order is character- 
ized by broken spin rotational invariance (^(S'f) 2 ^ 7^ 0, 

a = x,y or z) accompanied by unbroken time reversal 
symmetry ((Si) = 0). 

The most dramatic effect of the introduction of disor- 



der is seen in the BG phase, see Fig. 12(b). In that case, 
the presence of sites with different average occupations, 
both even and odd (see Fig. 9(a)), leads to the settling of 
the spatial average of Sf at a value intermediate between 
and 2. There is actually a very smooth dependence of 
this spatial average on t as we move from the SF into 
the BG phase. In the MI phase, on the other hand, the 
spatial average of Sf still vanishes and in the SF it also 
retains its generic intermediate values. 

The probability distribution of the expectation value 
of the square of the total spin per site Ps^(Sf) is shown 
in Fig. 13 for the two values of (x/Uq = 0.1 and 1.0 and 
several values of the hopping amplitude. From the pre- 
vious discussion, the behavior of Ps^(Sf) is expected to 
track closely the distribution of the average site occupa- 
tion P n (rii). Indeed, upon approaching the BG from the 
SF as t decreases, as shown in Fig. 13(a), the distribution 
of Sf becomes increasing broader with a bimodal shape, 
indicating the gradual appearance of both spin-zero and 
spin-1 sites, corresponding to the peaks at n,- = and 
rii = 1, respectively, of Fig. 9(a). Likewise, as t decreases 
and the system transitions from the SF to the MI, as 
depicted in Fig. 13(b), the Sf distribution shifts weight 
from non-zero values spread around ~ 1 down to a delta 
function at zero, at the same time as the average occu- 
pation distribution narrows down to a delta function at 
occupation = 2. 

It was argued in reference 32 that the spatially- 
averaged value of (Sf) intermediate between and 2 of 
the BG phase shown in Fig. 12(b) is indicative of a spin 
nematic phase 28 . In the clean case discussed in 28, how- 
ever, it is quantum inter-site spin correlations that are 
responsible for the appearance of nematic order. This oc- 
curs even at perturbatively small t, in which case each site 
has a fixed odd number of bosons. On the other hand, no 
inter-site spin correlations are incorporated in the SMFT 
and the BG phase is characterized by the presence of 
sites with different number occupations, see Fig. 9(a). It 
is these sites, with an odd number of bosons and Si = 1 
or an even number and Si — 0, which are ultimately re- 
sponsible for the intermediate value of (Sf) which both 
the SMFT and the Gutzwiller approach of reference 32 
find. This situation is rather different from the clean 
nematic and should not, in our view, be confused with 
it. 



F. The quantum phase transition as a function of 
disorder 

In previous Sections we fixed the disordered strength 
and analyzed the behavior of the system as a function of 
the chemical potential and the hopping amplitude. For 
the value of disorder we used (A/Uq = 0.3), the clean 
even-numbered MI lobes survived the introduction of 
randomness, whereas the odd-numbered ones were com- 
pletely destroyed. It would be interesting to see how 
the former behave as the disorder strength is further in- 
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Figure 13. (Color online) The probability distribution func- 
tions P S 2 (5f ) of the square of the total spin per site for var- 
ious values of the hopping amplitude and two values of the 
chemical potential: (a) n/Uo = 0.1 and (b) fi/Uo = 1.0. The 
disorder is set to A/Uo = 0.3. 



Figure 14. (Color online) Probability distributions functions 
of: (a) the order parameter and (b) the standard deviation 
of the occupation, for various values of disorder parameter 
A/Uo. The chemical potential is fixed at /j,/Uo — 1 and the 
hopping amplitude at t/Uo = 0.075. 



creased. We take up this task in this Section. 

We show in Fig. 14 the various distribution functions in 
the strong disorder regime (A > Uq) for h/Uq = 1.0 and 
t/Uo = 0.075. We remind the reader that for A = 0, this 
corresponds to a point well inside the rij = 2, Si = MI 
lobe. The order parameter distribution Pi (ipi) is shown 
in Fig. 14(a). For A/U Q = 1.0 the MI lobe has been sup- 
pressed in favor of the disordered SF phase, characterized 
by finite SF order parameters. As A is increased, this 
distribution broadens with increasing weight at ipi = 0. 
Eventually, for large enough randomness, the distribu- 
tion collapses to a delta function at ipi — 0, signaling 
the destruction of the SF phase. The distribution of 
site occupation standard deviation Pau (An^), depicted 
in Fig. 14(b), shows a qualitatively similar behavior, as 
expected. As discussed in Section IV D, this quantity 
largely tracks the distribution of the order parameters. 
But is the non-SF phase at large values of A a BG or a 
MI? 

One possible diagnostics tool is afforded by the dis- 
tribution of the mean square of the total spin per site, 
as shown in Fig. 15(a). It shows the typical broad, bi- 



modal shape characteristic of the BG distribution as A 
increases, indicating the presence of both singlet (weight 
at 0) and spin-1 (weight around 2) composites at each 
site. This should be compared with the similar small 
t distributions of Fig. 13(a), characteristic of the BG 
phase, and contrasted with the corresponding curves of 
Fig. 13(b), which are associated with MI behavior. 

Even more significant is the behavior of the distribu- 
tion of the average site occupation number P n (rij), shown 
in Fig. 15(b). As A increases, it is clearly seen that 
the distribution gradually evolves into essentially isolated 
peaks centered around the integer values (1 through 4 for 
A/Uq = 16). The presence of sites with different integer 
occupations in a non-SF phase is the hallmark of the BG, 
cf. Fig. 9(a). 

In order to dissipate any doubt that the large disorder 
phase for this particular choice of parameters is indeed a 
BG, we show in Fig. 16 the compressibility as a function 
of disorder strength. Its value decreases with increasing 
disorder but remains finite at the largest value analyzed 
(A/Uq — 16) at which point the order parameter has 
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Figure 15. (Color online) Probability distributions functions 
of: (a) the expectation value of the square of the spin op- 
erator per site and (b) the average site occupation number, 
for various values of disorder parameter A/Uq. The chemical 
potential is fixed at fi/Uo = 1 and the hopping amplitude at 
t/U = 0.075. 



already vanished, cf. Fig. 14(a). 

This shows that, though more robust against random- 
ness, the even-numbered MI lobes can also be wiped out 
and transformed into BG phases with sufficiently large 
disorder. It follows that, for large enough values of A, 
only the SF and the BG phases survive, as had been 
previously observed in the spin-zero case ' 35 . 



V. CONCLUSIONS 

We have analyzed the stochastic mean field theory of 
the disordered spin-1 Bose-Hubbard and discussed its 
physical properties as a function of the hopping ampli- 
tude, the chemical potential and the disorder strength. 
Although the model exhibits strong similarities with 
its spin-zero counterpart, several differences stand out. 
There is a clear difference in the behavior of the odd- 
and the even- numbered MI lobes. The latter are much 
more robust with respect to the introduction of disorder. 
As a result, there is a sizable portion of the parameter 
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Figure 16. (Color online) The compressibility as a function of 
disorder strength. The chemical potential is fixed at (j,/Uo = 1 
and the hopping amplitude at t/Uo = 0.075. 



space in which only even-numbered MI lobes exist, the 
odd- numbered ones having been transformed into a BG. 
The BG insulator is characterized by a finite compress- 
ibility and the presence of sites with different occupa- 
tions, like the spin-zero case. However, unlike the spin- 
zero case, different occupations give rise to different spins. 
Therefore, the spin-1 BG is an inhomogeneous mixture 
of spin-0 and spin-1 composites within the SMFT. Very 
similar behavior was obtained within the Gutzwiller ap- 
proach of reference 32. We should stress that reference 
32 employs two different approaches in the study of the 
disordered system. In one approach, which the authors 
call a 'probabilistic mean-field theory', only an average 
order parameter is considered in the description. This 
approach is a much poorer description than the present 
SMFT, since it incorporates no spatial fluctuations and, 
in particular, does not exhibit a BG phase. Alterna- 
tively, they also show a direct lattice calculation of the 
site-decoupled mean-field theory. This does have spatial 
fluctuations and describes all 3 phases. It includes spa- 
tial correlations of the local order parameter which are 
absent in our SMFT treatment and should therefore be 
considered a superior approach. However, direct com- 
parison shows that the phase diagram and some physical 
properties we obtain are almost exactly the same as the 
lattice calculation of reference 32, highlighting that the 
much simpler SMFT already incorporates the most im- 
portant correlations of the more complete treatment. 

The presence of local composites with different total 
spin values raises the important question of the spin cor- 
relations within the spin-1 BG phase. As discussed in 
28 for the clean system, spin correlations outside the 
scope of either the 'site-decoupled' mean-field theory or 
the Gutzwiller approach give rise to non-trivial nematic 
ground states. It would be of great interest to incorporate 
these into a description of the disordered system in or- 
der to investigate the interplay between disorder-induced 
number fluctuations and quantum inter-site correlations. 
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Besides the possibility of a spin nematic, the introduc- 
tion of randomness could also potentially give rise to 
spin-glass order, a Bose-spin- glass or quantum Griffiths 
phases 45-48 . 

Another direction deserving of further scrutiny is the 
case of ferromagnetic interactions, U2 < 0. In this case, 
we expect the MI lobes to be characterized by the bosons 
aligning to form a maximum spin composite. In the pres- 
ence of strong enough disorder, spins of different sizes are 
expected to form, rendering the problem of the ground 



state spin structure even richer. 
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